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The  stability  theory  for  finite  difference  Initial  Boundary- Value  approximations  to  sys¬ 
tems  of  hyperbolic  partial  differential  equations  states  that  the  exclusion  oLeigenvalues  r.nd 
generalized  eigenvalues  is  a  sufficient  condition  for  stability.  The  theory,  however,  does  not 
discuss  the  nature  of  numerical  approximations  in  the  presence  of  such  eigenvalues. 

In  fact,  as  was  shown  previously^,  for  the  problem  of  vortex  shedding  by  a  cylinder 
in  subsonic  flow,  stating  boundary  conditions  in  terms  of  the  primitive  (non-characteristic) 
variables  may  lead  to  such  eigenvalues,  causing  perturbations  that  decay  slowly  in  space  and 
remain  periodic  time.  Characteristic  formulation  of  the  boundary  conditions  avoided  this 
problem. 

In  this  paper,  we  report  on  a  more  systematic  study  of  the  behavior  of  the  (linearized) JJ 
one-dimensional  gas  dynamic  equations  under  various  sets  of  oscillation-inducing  [‘legal"  ' 
boundary  conditions.  /  /  ,  j 

)  t —  v 
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1.  Introduction 


The  increase  in  computers’  speed  and  memory  allows  researchers  to  investigate  fluid 
dynamical  problems  with  greater  attention  to  delicate  features  of  the  flows.  For  example, 
only  in  the  last  five  years  have  investigators  [1,6]  computed  the  vortex  shedding  phenomenon 
behind  a  two-dimensional  cylinder  at  low  Mach  numbers  and  moderate  Reynolds  numbers. 
Unfortunately,  with  this  improved  ability  due  to  increased  computer  power  to  compute 
complex  phenomena  comes  the  necessity  to  deal  with  unforeseen  numerical  surprises  which 
might  be  mistaken  for  real  or  physical  effects.  For  example,  in  the  case  of  the  2-D  cylindrical 
Von-Karman  vortex  street  computation  -  the  main  features  of  the  flow,  such  as  the  shedding 
frequency  are  predicted  accurately.  However,  a  concomitant  computational  result  appearing 
in  the  output  data  base  is  a  spurious  secondary  frequency  which  was  mistakingly  attributed 
to  the  start  of  transition.  In  previous  work,  [1],  we  have  shown  that  the  spurious  secondary 
frequency  resulted  from  applying  the  far-field  inflow  boundary  condition  to  the  primitive 
variables.  It  was  also  shown  there  that  the  same  boundary  treatment,  applied,  however, 
to  the  characteristic  variables,  eliminated  this  phenomenon.  In  this  paper,  we  characterize, 
by  an  analytic  description,  the  boundary  conditions  (both  inflow  and  outflow)  under  which 
the  numerical  solution  of  the  Euler  equations  will  exhibit  temporal  oscillations  which  are 
foreign  to  the  exact  solution  of  the  p.d.e.’s.  The  starting  point  of  our  study  is  the  modal 
analysis  developed  by  G-K-S  [2],  and  Osher  [3].  This  theory  states  that  for  finite  differences 
approximations  to  Initial  Boundary  Value  problems  of  hyperbolic  systems  stability  is  assured 
by  the  exclusion  of  eigenvalues  and  generalized  eigenvalues.  The  theory,  however,  does  not 
discuss  the  nature  of  the  numerical  approximations  in  the  presence  of  such  eigenvalues. 

The  model  that  we  study  is  that  of  the  1-D  compressible  Euler  equations  linearized  about 
free  stream  conditions.  The  numerical  schemes  that  we  analyze  are  second  order  in  space 
and  time  represented  by  the  Lax-Wendroff  scheme.  In  our  1-D  case  it  is  equivalent  to  all 
other  second  order  algorithms  such  as  the  MacCormack  scheme  [4], 

In  Section  2,  we  briefly  review  the  underlying  theoretical  considerations  both  for  the 
p.d.e.  and  the  f.d.e.  formulation. 

In  Section  3,  we  discuss  the  case  of  inflow  b.c.’s  expressed  with  primitive  variables.  We 
show  that  even  though  the  numerical  solution  is  technically  stable  and  therefore,  by  the  Lax 
equivalence  theorem,  convergent,  for  finite  meshes  and  time  it  shows  spurious  oscillations 
that  decay  only  slowly  with  mesh  refinement.  In  Section  4,  we  repeat  this  demonstration 
for  the  case  of  primitive- variable  description  of  the  out-flow  boundary  conditions.  In  Section 
5,  we  apply  a  primitive-variable  formulation  of  the  boundary  conditions  both  at  the  inflow 
and  outflow  boundaries.  This  allows  a  nonlinear  interaction  between  the  modes  created  at 
both  boundaries.  At  this  point,  the  reader  should  be  reminded  that  the  same  treatment 
of  the  boundary  conditions,  applied  to  the  characteristic  variables  eliminates  any  spurious 
frequencies.  This  is  a  corollary  of  our  analytical  formulation  and  is  borne  out  in  our  numerical 
experiments. 

All  sections  contain  numerical  examples  pertinent  to  the  initial  boundary  value  problem 
considered  therein. 
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2.  Analytical  Preliminaries 


We  start  with  the  one-dimensional  Euler  equations  of  gas  dynamics  in  conservation  form: 

dU  dFVU)  n  ,  % 

1  A  '  ~ n  (2.1) 


dt  dx 


=  0 


where 


U  = 


F  = 


P 

m 

E 

m 


771  i  __ 

—+p 


(2.2) 


(2.3) 


L  u(E  +  p)  J 

Here  u,  p,  p  and  E  are,  respectively,  the  velocity,  density,  pressure,  and  total  energy  per  unit 
volume,  m  =  pu  is  the  mass  flux.  In  the  case  of  an  ideal  gas  the  equation  of  state  is 


P=  (7  ~  1) 


„  1 

(  m2\l 

E  -  - 

2 

V  P  )\ 

(2.4) 


where  7  =  7/5  for  diatomic  gases  such  as  air.  Equation  (2.1)  can  be  rewritten  in  non¬ 
conservation  form  as 

f  +  =  °  <2-5) 

where  A  =  dF/dU  is  the  Jacobian  of  the  flux  vector  F  with  respect  to  the  solution  vector 
U.  Linearizing  about  steady  free  stream  conditions  ,  U £  =  (poo,  PooWoo,  -£«>),  Equation  (2.5) 
becomes: 

gj(«0  +  M.U^{SU)  =  0,  (2.6) 

where  6U  =  U  —  {/<*,  is  the  perturbation  vector.  The  matrix  A{Uoo )  =  has 

three  eigenvalues,  aj  —  —  c^,  a2  =  u ^  +  cx  and  a3  =  u^.  The  free  stream  sound 

speed,  Coo,  is  given  by  Coo  =  (iPoo/Pao)1-  The  corresponding  eigenvectors  in  terms  of  the 
conserved-variable  perturbations  are  given  by 


Ri  =  -{6m  -  UooSp)  + 
R2  —  {6m  —  Uoo 6p)  + 


7  ~  1 

Cqq 

7  - 1 


[\uloSP  -  uoo 6m  +  6E 

\ulo^P  ~  Woo 6m  +  6E 

,  £ 


and 


R3  -  Co,  6 p 


7  -  1 


\ulo^P  ~  Woo 6m  +  6E 


Furthermore,  using  a  linearized  version  of  the  equation  of  state,  i.e., 

n 


6p  =  {l~  1) 


^u^p  -u^mF  6E 


(2.7) 

(2.8) 

(2.9) 

(2.10) 
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equations  (2.7)  -  (2.9)  may  also  be  written  in  terms  of  the  primitive- variable  perturbations, 
( 6p,6u,Sp ): 

Ri  =  -PooSu-\ - 6p,  (2-11) 


Ri  =  PooSu  -j - Sp, 


(2.12) 


R3  =  - <5p. 


(2.13) 


Any  perturbation  imposed  on  the  free  stream  solution  will  evolve  as  a  combination  of  these 
eigenvectors. 

In  terms  of  the  characteristic  variables  Rt(a  =  1,2,3)  defined  in  (2.7  -  2.9)  and/or  (2.11 
-  2.13),  Equation  (2.6)  may  be  written  as: 


dRi  ,  ,  sdR!  _  n 

Qt  +  (ti»  Co.)  gx  ~  0 

dR2  f  .  dR2  _ _ 

dt  +  ^°° +  c°°^  dx  ~ 0 
dR3  dR3 

-ar+“”-§r  =  0 

j.  „  dR-  -  a  \  i 

at  +  ‘ax  ~0,  1’2,3’ 


(2.14) 


(2.15) 

(2.16) 


(2.17) 


In  a  finite  domain,  say  0  <  x  <  1,  for  the  subsonic  case  <  c^,  the  system  (2.7)  is  well 
posed  with  the  following  initial  and  boundary  conditions: 


Ra(x,0)  =  fa(x)  3  =  1,2,3 


(2.18) 


R?{0,t)  -  aoRi{0,t)  =  9l{t)  (2.19) 

■^3(0,  t)  —  0oRx(O,  t)  =  <72(2)  (2.20) 

•Ki(l>£)  4-  01  #2(1, 0  +  £1^3(1, t)  =  53(0  (2  21) 

where  ac,,f30,ai  and  E\  are  arbitrary  real  constants  and  f,{x),g,(t )  (s  =  1,2,3)  are  square 
integrable  in  their  respective  domains. 

We  get  numerical  approximations  of  second  order  spatial  and  temporal  accuracy  by  using 
the  Lax-Wendroff  scheme, 

j>0  ,  =  1,2,3  (2.22) 

where  w*,n  =  w‘(j Ax,nAt)  is  the  finite  difference  approximation  1.0  R,(x,i). 

We  note  at  this  junction  that  although  we  shall  illustrate  tb^  detailed  development  of  the 
spurious  frequencies  using  the  Lax-Wendroff  scheme,  they  depend  in  fact  mostly  on  the  form 
of  the  finite-difference  boundary  conditions  (yet  to  be  specified)  and  not  on  the  particular 
inner  algorithm. 
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To  close  the  system  (2.22)  we  need  both  the  finite  difference  form  of  the  algebraic  bound¬ 


ary  conditions  (2.19  -  2.22)  i.e., 

w20'n  -  aawo'n  =  ?i(nAt)  (2.23) 

w30'n  -  6jWl0'n  =  g2(nAt)  (2.24) 

w#n  x  +  £iw3tin  =  g3(nAt)  (2.25) 

and  also  three  “numerical”  boundary  conditions,  one  at  x  =  0  and  two  at  x  =  1  =  A xN.  A 
commonly  used  method  of  imposing  “numerical”  boundary  condition  is  to  extrapolate  from 
the  interior  in  the  following  manner: 

w0'  +  gqWq  -f  CqWq  —  w{  +  a0w{  -|-  £0^i  >  (2.26) 

and  at  the  outflow  boundary,  x  =  1  =  N Ax, 

w2jf  -  aiw]}n  =  w2f}\  -  a (2.27) 
wn  ~  Piw]}n  =  wn- i  —  (2.28) 


The  “numerical”  boundary  conditions  (2.26)  -  (2.28)  are  zeroth  order  extrapolations,  with 
ai,(3i,a0  and  e0  arbitrary  real  constants. 

For  stability  analysis  studies  it  is  sufficient  to  consider  the  case  of  homogeneous  boundary 
conditions,  i.e.,  gT(nAt)  =  0 ,r  =  1,2,3.  We  look  for  solutions  of  the  form 

w‘'n  =  zn(A,K\  +  s  =  1,2,3.  (2.29) 

Substituting  this  ansatz  into  (2.22)  we  find  that  the  k,'s  and  fi,’ s  are  the  roots  of  the 
quadratic  equation 

(A,2  -  \.)x)  +2(1  -  *  -  \])X.  +  (A2  +  A.)  =  0  (2.30) 


where  X,  =  a, At/ Ax.  For  \z\  >  1  one  of  the  roots  of  (2.30)  (say  k,)  is  inside  the  unit  circle 
and  the  other,  g.,,  is  outside  the  unit  circle. 

Substituting  (2.29)  into  the  homogeneous  version  of  (2.23)  -  (2.28)  we  get  the  following 
system  of  equations: 

/  A:  \ 

a2 

Q(k,(z),h,(z))  g3  =0  (2.31) 

B2 

\  B3  ) 
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where  Q(k,(z),  p,(z))  =  Q(z)  is  given  by 


/Cl  —  1 

cr0(K2  -  1) 

£o(«3  -  1) 

Mi  ~  1 

Oo(M2  -  1) 

£o(M3  -  1) 

-0-0 

1 

0 

-a0 

1 

0 

—00 

0 

1 

-0o 

0 

1 

kn 

K1 

01*2 

£lK3 

„N 

Mi 

£lM3 

-<*l>cf-1(Kl  - 

1) 

«f_1(K2  -  1) 

0 

-o=iMf-1(Mi 

-1) 

M^~1(M2  -  1) 

0 

.  -01K?~1(/<1  - 

1) 

0 

-  1) 

-1) 

0 

M^'Hms  -  1) 

(2.32) 

Note  that  the  dependency  on  z  in  (2.32)  comes  from  the  fact  that  (j  =  1,2,3)  are 

functions  of  z  defined  implicitly  in  (2.30). 

It  is  quite  easy  to  formulate  a  necessary  condition  for  stability.  This  is  the  Riabenkii- 
Godunov  condition  [5]: 

Lemma  2.1:  The  Lax-Wendroff  scheme  (2.22)  with  the  inflow  boundary  conditions  (2.23), 
(2.24),  and  (2.26)  and  outflow  conditions  (2.25),  (2.27),  (2.28)  is  unstable  if  there  is  \zq\  >  1 
such  that 

det  Q(z0 )  =  0.  (2.33) 

In  fact,  such  a  zo  gives  a  solution  of  the  type 

=  *o{A.k{  +  B.n>) 

that  grows  with  the  number  of  time  steps.  This  is  indeed  a  case  of  classical  instability. 

The  sufficient  condition  for  stability  can  also  be  formulated  in  terms  of  the  determinant 
of  Q(z );  in  fact  it  has  been  proven: 

Lemma  (2.2):  [GKS]  The  LW  scheme  (2.22)  with  the  b.c.’s  (2.23)  -  (2.28)  is  stable  if  for 
every  \z0\  >  1, 

det  Q{z0)  ±  0  (2.34) 

It  is  evident  that  the  case  |20|  =  1  with  det  Q(z0)  =  0  is  covered  by  neither  lemma.  It  is 
precisely  this  case,  however,  that  is  responsible  for  the  spurious  oscillations. 

For  the  sake  of  convenience,  we  rewrite  (2.31)  and  (2.32)  in  the  following  way: 
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Define 


B(k) 


Now  (2.31)  becomes 


and 


A(k)  = 


I  /Cl  -  1  Co (k2  -  1)  £0(«3  -  1J  \ 
1  0 


— a0 


V  -Po 


0 


1 


/ 


K 


—  OCiKi  1(«1  —  1)  2  —  1) 


N 


n  kN 
<7l  K2 


N~l, 


F  KN 


^  1(K1  -  1) 


A(k)  a  +A(fi)  6=0 
B(ic)  a  +B(fj,)  6=0 


Q(z) 


A(k)  A(h)  \ 
B(k)  B(h)  ) 


(2.35) 


(2.36a) 

(2.366) 

(2.37) 


The  determinant  condition  (2.33)  involves  the  solution  of  a  very  complicated  nonlinear  com¬ 
plex  equation.  Only  little  insight  can  be  gained  by  trying  to  directly  analyze  it.  One  simple 
case,  though,  is  given  by  the  fully  characteristic  case 


c*o  —  00  =  a0  =  e0  —  a1  ~  f3i  —  0\  —  C\  —  0. 
Note  that  in  this  case  A  and  B  are  diagonal 


/  /cx  -  1  0  0  ^ 

o 

o 

A(k)  = 

0  1  0 
i  o  0  lj 

B(k)  = 

0  k2  1(k2  -  1)  0 

^o  0  k%-\k3-1  )} 

(2.38) 


(2.39) 


and  the  determinant  condition  reduces  to 


[(«i-l)Mf-«f(Mi-l)][(Ma-l)^  1(«2  — l)][(/x3  — l)/x^  1(«3-1)]  =  0.  (2.40) 

Equation  (2.40)  indicates  that  the  choice  of  parameters  in  (2.38)  decouples  the  system  (2.31) 
and  the  system  is  thus  reduced  to  the  scalar  case.  It  may  be  easily  shown,  using  (2.30),  that 
there  is  no  zq  for  which  (2.40)  has  a  solution,  and  therefore  no  spurious  solutions  exist  in 
the  fully  characteristic  case. 
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3.  The  Inflow  Case 


In  this  section,  we  will  investigate  the  possibility  of  having  a  spurious  frequency  induced 
by  “primitive”  boundary  conditions  (numerical  and  analytical)  at  the  inflow  boundary.  By 
“primitive”  boundary  conditions,  we  mean  b.c.’s  expressed  in  terms  of  primitive  variables. 
By  a  spurious  frequency,  we  mean  a  solution  of  the  type  z0  =  el<t>  to  the  determinant  condition 
(2.33).  Indeed,  if  the  determinant  has  such  a  root  zo,  then  a  spurious  solution  of  the  form 

W’n  =  zn(A.Kl  +  B.ni)  =  e'^A.4  +  B.p?,)  (3.1) 

is  superimposed  on  any  numerical  solution  of  (2.22)  -  (2.28). 

When  the  number  of  grid  point  increases  without  limit,  N  — »  oo,  the  influence  of  the 
outflow  boundary  treatment  gets  decoupled  from  that  of  the  inflow  boundary.  This  can  be 
seen  by  observing  that  the  submatrix  B(/c)  defined  in  (2.35)  vanishes  because  |/c|  <  1.  At 

the  same  time  I?(/x)  tends  to  infinity.  Therefore,  from  (2.36b)  we  have  that  b  — ►  0.  We  have 
then,  from  (2.36a) 

A(/c)o=0  (3.2) 

and  the  determinant  condition  (2.32)  reduces  by  (2.37)  to 

det  Q(z )  =  (*x  -  1)  +  a0a0(K2  -  1)  +  £o/3o(«3  -  1)  =  0.  (3.3) 

In  practice,  one  cannot  achieve  N  — »  oo.  We  then  inquire  how  to  best  decouple  the 
outflow  influence  from  the  inflow  for  the  case  of  fixed  N,  albeit  large.  If  we  set  the  amplitude 
of  the  perturbations  due  to  outflow  treatment,  |  b  |,  equal  to  zero  then  from  (2.36)  we  have 

(A(/c)  +  B(k))  o=  0  (2.36c) 

as  well  as  A(k)  a=  0.  Thus  the  term  B(k)  in  (2.36c)  represents  the  effect  of  the  incomplete 
decoupling  due  to  the  finite  N.  Since  N  1,  it  is  clear  from  the  definitions  that  B(k) 
represents  a  small  perturbation  to  A{k),  which  we  would  like  to  minimize.  It  can  be  shown 
that  this  minimization  is  achieved  for  £j  =  crj  =  ai  =  /3\  =  0.  Notice  that  this  corresponds 
to  imposing  both  analytic  and  numerical  outflow  boundary  conditions  on  the  characteristic 
variables,  Rt(3  =  1,2,3). 

A  spurious  solution  exists  if  the  determinant  condition  (3.3)  has  a  root  of  the  form 
z0  =  e'*.  In  this  case,  the  homogeneous  version  of  (2.21)  -  (2.28)  with  =  ex  =  aj  =  fli  —  0 
admits  a  solution  of  the  form 

w’j'n  =  A,exn+i cJ„  <f>  ±  0 

<  1- 

Several  observations  can  be  made  at  this  stage: 

(i)  No  spurious  frequency  is  created  if  either  the  analytical  inflow  boundary  conditions 
(2.23)  -  (2.24)  ot  the  numerical  one  (2.25)  are  in  characteristic  form. 

If  the  analytical  boundary  conditions  are  in  characteristic  form,  then  a0  —  /30  =  0  and 
the  system  case  is  reduced  to  the  scalar  case.  If  the  numerical  inflow  condition  (2.26) 
is  in  characteristic  form,  then  cr0  =  e0  =  0  and  the  system  again  reduces  to  the  scalar 


(3.4) 

(3.5) 


(ii)  On  the  other  hand,  if  the  inflow  boundary  conditions  are  not  given  in  terms  of  char¬ 
acteristic  variables,  then  there  is  a  wealth  of  possible  boundary  treatments  that  lead 
each  to  different  spurious  solutions.  In  fact,  for  every  frequency  there  exists  a  set  of 
inflow  boundary  conditions  that  induces  a  spurious  solution  with  that  frequency.  In 
particular,  given  an  arbitrary  Zq  =  e’*;  <f>  ^  0  we  get  /c,(zo)  from  (2.30).  We  can  then 
view  (3.3)  as  an  equation  for  ao <?o  an<^  A)£o-  Since  these  quantities  are  real,  whereas 
the  k,(z0)  are  complex,  Equation  (3.3)  is  a  system  of  two  equations  for  Q0a0  and  S0e0, 
yielding  a  unique  solution.  The  explicit  form  of  the  spurious  solution  (3.4)  reveals  the 
nature  of  such  a  solution.  We  summarize  the  results  in  the  results  in  the  following 
lemma: 


Lemma  (3.1):  The  spurious  solution  ( S.f )  converges  to  zero  for  every  fixed  x  =  jAx  and 
t  =  n/\t  as  Ax  — ♦  0,  At  — >  0. 


Proof:  Since  Ax  — >  0  and  x  is  fixed,  then  j  — »  oo.  Note  that  |k,|  <  1,  and  therefore 
lim.,_o  =  0. 


Notice  that  even  though  Lemma  (3.1)  indicates  convergence  of  the  total  scheme,  still  on 
a  fixed  finite  grid,  the  spurious  solution 


U); 


may  be  confused  with  a  real  time  periodic  solution.  In  particular,  if  one  of  the  /c,(e"*’)  is 
close  to  unity  in  magnitude,  then  the  time  periodic  solution  may  show  up  in  a  large  part 
of  the  spatial  domain.  Of  course,  if  all  the  k’s  are  small  in  magnitude,  the  time  periodic 
spurious  solution  will  be  confined  to  a  very  narrow  boundary  layer  at  the  inflow  boundary. 


To  illustrate  the  above  analysis,  we  chose  to  solve  numerically  equations  (2.14)  -  (2.16) 
with  Uoo  —  1  and  M «,  =  .4.  Thus,  ai  =  — 1.5;o2  =  3.5;  a3  =  1 ;  c^  =  2.5.  For  the  initial 
conditions  (2.18),  we  chose 

J2i(*,0)  =  e"1  J2a(x,  0)  =  eI  i23(*,0)  =  e2x.  (3.7) 

For  the  analytical  boundary  conditions  (2.19)  -  (2.21),  v*e  chose 

a0  =  .582155;  p0  =  -.6115;  cti  =  0;  ex  =  0  (3.8) 

and 

Si(0  =  e“°2t  -  a0eait 

g2(t)  =  e-2a’t-/VM  (3.9) 

ffa(f)  =  e°lt. 

It  is  readily  verified  that  under  the  above  conditions  the  exact  solution  for  the  system  is 

Ri(x,t)  =  e-(x~ait) 

R2(x, t)  =  ex~a^  (3.10) 

#3(x,0  =  e2(*-°st). 
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Thus  the  analytic  solution  decays  to  zero  at  steady  state.  The  system  was  discretized  by 
the  Lax-Wendroff  scheme  (2.22)  with  the  following  parameters  for  the  numerical  boundary 
conditions  (2.26)  -  (2.27)  -  (2.28) 


=  A  =  0 

cq  —  1  Cq  =  2. 


(3.11) 


Note  that  at  the  outflow  boundary,  x  =  1,  both  the  analytic  and  the  numerical  boundary 
conditions  are  in  characteristic  form,  because  in  this  section  we  are  interested  in  the  influence 
of  inflow  b.c.’s  only. 

It  can  be  shown  that  the  inflow  determinant  condition  (3.3)  has  for  the  above  choice  of 
parameters  a  solution  of  the  form 

z0  =  e°'3i.  (3.12) 

The  corresponding  /cf(zo)  obtained  by  the  quadratic  equation  (2.30)  satisfy 


|/ci|  =  . 23035  |/c2|  =  . 95334  |/c3|  =  .19015  (3.13) 


In  the  following  graphs,  we  present  w”’s  for  several  meshes 


A  -  1  1  1  1 

X  _  16’ 32  ’  64’ 128 


(3.14) 


At  was  taken  to  be  O.lAx,  corresponding  to  a  CFL  of  0.35  for  our  flow  parameters.  The 
results  are  sampled  at  x  =  .5,  0<<<9;  thus 


j Ax  =  .5  and  0  <  nAt  <  9.  (3.15) 

At  t  =  9  the  analytic  solutions  R,  (|,9)  (s  =  1,2,3)  are  of  the  order  of  10~6. 

The  numerical  solution  uyn  versus  time  is  presented  for  x  —  .5  and  several  meshes  in 
figures:  l(a,b,c,d)  for  A-  =  N  =  16, 32, 64, 128.  We  do  not  present  and  w3,n  because  they 
converge  even  on  coarse  meshes  ( N  >  16)  to  the  exact  solution.  This  convergent  behavior  at 
x  =  .5  is  due  to  the  fact  that  |/ci|  and  ]/c3|  are  small  and  hence  |«i|^2  and  |/c3 are  already 
of  the  same  order  (~  10-6)  as  the  analytic  solution.  On  the  other  hand,  displays  a 
time  periodic  solution  on  each  of  the  meshes  presented.  However,  as  expected,  as  the  mesh 
is  refined  the  spurious  period  is  halved  and  the  amplitude  decreases.  The  solution  converges 
with  the  decreasing  mesh,  but  for  every  one  of  the  fixed  meshes  used  the  solution  at  x  =  .5 
still  exhibits  the  spurious  frequency,  because  \k2\n^2  is  not  yet  small  enough.  This  behavior 
at  x  =  .5  is  representative  of  the  solution  at  other  spatial  locations. 


4.  The  Outflow  Case 


In  this  section,  we  deal  with  the  effects  of  non-characteristic  formulation  of  the  outflow 
boundary  conditions  in  a  similar  manner  to  that  of  Section  3. 

It  is  easier  to  do  the  asymptotic  ( N  — »  oo)  analysis  in  the  outflow  case  by  using  new 
coefficients  C’s.  Defining 


c  — 


(  Cl\ 
C2 

\  ^3  ) 


(4.1) 
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we  may  rewrite  (2.36)  as  follow-1- 


A(k)S 4-  —  0 

B(k)S  +  E(fi)c  =  0 


(4.2) 

(4.3) 


where  j 

D(im)  =  A(/x)A(/x)  (4.4)  ; 

E{n)  =  B(fi)A(ii)  (4.5)  : 

(4-6) 

/  I 

and  the  matrice  A  and  B  are  defined  in  (2.35). 

Again,  for  N  large,  D([i)  and  £(/x)  have  small  entries.  In  order  to  decouple  the  influence 
of  the  inflow  conditions,  we  shall  assume  a=  0  and  so  the  determinant  condition  becomes 


A  in)  = 


Mi 


N+1  0  0 
0  ^ N+1  0 
0  0  ^N+1 


or,  more  explicitly 


det  j£(/i)  =  0 


(4.7) 


det  Ein)  =  Ati(/x2  -  1)(^3  -  1)  +  a1cr1fta(/x1  -  l)(/i3  -  1)  + -  1)(M2  -  1)  =  0.  (4.8) 
A  technical  argument  similar  to  the  one  leading  to  (3.4),  (3.5)  gives  a  solution  of  the  form 


w’-n  =  C,e'n“ix?-i,  uj^O 

\n.(e'un)\  >  1 


(4.8) 


The  observation  made  there  are  also  valid  here,  namely  that  a  spurious  solution  (as  defined 
there)  may  exist  for  non-characteristic  outflow  treatment;  see  observations  (1)  and  (2)  on 
pages  7-8.  The  numerical  examples  are  again  based  on  (2.14)  -  (2.16)  with  the  same  flow 
parameters,  and  with  the  same  initial  conditions  given  by  (3.7).  For  the  analytic  boundary 
conditions,  we  took  a0  =  /30  =  0,  a i  =  l,£i  =  2,  and,  from  (2.23)  -  (2.25), 
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92 


e~°st 

g  — 2aj( 

eait  +  <r_  e~a,t  +  eie~2a3t. 


The  analytic  solution  is  again 

Riix,t)  =  e-(x-ai‘),  R3ix,t)  1  =  ex-°3t,  R3(x,t)  =  e2(x-°st).  (4.9) 

The  numerical  boundary  conditions  (2.26)  -  (2.28)  imposed  on  (2.22)  were: 


at  =  -4.668,  fit  =  2.09485,  e0  =  a0  =  0.  (4.10) 

These  values  of  c*i  and  fit  were  chosen  so  that  we  will  have  a  case  of  a  spurious  freauency 
accompanied  by  at  least  one  /x  near  to  1.  The  resulting  /x’s  and  z0  are: 

|/ii  j  =  1.00625,  |/i2|  =  2.0774,  |/i3|  =  1.3959,  z0--eou.  (4.11) 


The  numerical  results  for  w are  presented  (again  at  x  =  .5)  in  Figures  2(a,b,c,d,e).  w^J2 
and  display  a  convergent  behavior  due  to  the  smallness  of  \fi~  "+W,(s  =  2,3).  On 
the  other  hand,  io^"2  is  oscillatory  on  all  the  meshes  we  used,  although  the  amplitude  decays 
in  accordance  with  the  law. 

We  have  thus  demonstrated  that  spurious  frequencies  can  be  induced  by  improper  (i.e., 
non-characteristic)  treatment  of  boundary  conditions,  not  only  at  the  inflow  boundary  (see 
Ref.  1  for  the  physical  case  of  a  flow  past  a  cylinder)  but  also  at  the  outflow  boundary. 

5.  The  Finite  Domain  Case 

In  this  section,  we  report  on  numerical  experiments  in  which  both  boundaries  (x  =  0  and 
i  =  l)  are  treated  in  non-characteristic  fashion.  We  chose  the  same  boundary  conditions 
reported  on  in  the  previous  two  sections,  but  now  they  are  applied  simultaneously.  Recall 
that  when  only  the  inflow  was  non-characteristic,  the  solution  exhibited  spurious  oscillations 
but  the  solution  converged  as  Ax  — >  0,  At  — >  0.  The  outflow  case  on  its  own  exhibited  the 
same  behavior.  We  mention  also  that  the  G-K-S  theory  [2],  for  the  stable  case  (see  Lemma 

(2.2) ),  predicts  that  if  the  semi-infinite  cases  (i.e.,  pure  inflow  or  pure  outflow)  are  each 
stable  separately,  then  also  the  finite  domain  case  will  be  stable  and  convergent  with  mesh, 
though  perturbation  might  grow  in  time.  In  the  numerical  experimentation  reported  herein, 
it  will  be  seen  that  in  our  special  case  (|zo|  =  1  for  each  of  the  semi-infinite  problems)  the 
solutions  seem  to  grow  exponentially  in  time  and  to  converge,  but  non-uniformly,  with  mesh 
size. 

The  numerical  approximation  (2.22)  (for  s  =  1,2,3)  was  used  to  solve  (2.14)  -  (2.16) 
with  the  boundary  conditions  now  given  by: 

=  0.582155  A,  =  -0.6115  a0  =  1  e0  =  2  (5.1) 

=  -4.668  A  =  2.09485  ox  =  1  ex  =  2  (5.2) 

The  boundary  parameters  (4.1)  were  taken  from  the  “pure”  inflow  problem  and  those  in 

(4.2)  from  the  “pure”  outflow  problem. 

We  show  the  temporal  behavior  of  Wj'n  at  x  —  .5  in  Figures  3(a,b,c,d,e,f).  Examination 
of  these  figures  shows  that  unlike  the  case  of  “pure”  inflow  and  outflow,  where  there  were 
oscillations  but  they  did  not  grow  in  time  in  the  present  “finite  domain”  boundary  treatment 
the  oscillations  row  temporally.  w^'n  and  Wj,n  display  the  same  behavior  and  are  not  shown 
here.  Another  difference  between  the  “pure”  case  and  the  present  finite-domain  case  has 
to  do  with  behavior  as  Ax  is  being  decreased.  In  the  previous  sections,  we  saw  that  the 
amplitude  of  the  oscillation  decreased  uniformly  with  the  mesh  size.  In  the  present  case, 
there  is  an  increase  in  amplitude  as  Ax  is  decreased  from  N  =  16  to  N  =  128.  Afterwards, 
for  N  =  256  and  512  we  see  a  decreased  amplitude;  and  although  the  perturbing  oscillations 
are  still  large  compare  to  10-6  (and  growing  temporally)  there  does  seem  to  be  a  beginning 
of  convergence. 

We  should  remark  here  that  this  type  of  behavior  of  the  numerical  approximation  to  the 
solution  can  be  very  unsettling  for  practitioners  using  practical  codes,  since  multidimensional 
mesh  sizes  seem  likely  not  to  be  in  the  asymptotic  range  in  the  near  future. 
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